# ============================================================================
# Study 2 Replication Script
# Gender Differences in Immigration Attitudes in Japan 
# Japanese Sample (May 2024) - Starting from jp2_replication.RData
# ============================================================================

# Load necessary libraries
library(tidyverse)  
library(broom)      
library(ggeffects)  
library(scales)     
library(stargazer)  
library(forcats)    
library(gridExtra)  
library(knitr)      
library(kableExtra) 
library(texreg)     


# ============================================================================
# LOAD CLEANED DATA
# ============================================================================

# Load the cleaned dataset
load("study2_replication_data.RData")

# ============================================================================
# FIGURE 6: Predicted Values of Competition Concerns (MANUSCRIPT)
# In the published manuscript, these are presented as separate Figure 6a and 6b
# for clearer visualization of each outcome.
# ============================================================================

#Run Regression Models
# Model 1: Competition Concerns against immigrant women
model_women_concern <- lm(w_cncrn_r ~ conditions:female +
                            conditions + female +
                            inc_11_r +
                            edu_8_r +
                            age_5_r +
                            pride_r +
                            econ_nation_notsat_r + econ_self_notsat_r +
                            job_conf_7_r + 
                            pid_ldp,
                          data = jp2_replication) 

# Model 2: Competition Concerns against immigrant men
model_men_concern <- lm(m_cncrn_r ~ conditions:female +
                          conditions + female +
                          inc_11_r +
                          edu_8_r +
                          age_5_r +
                          pride_r +
                          econ_nation_notsat_r + econ_self_notsat_r +
                          job_conf_7_r + 
                          pid_ldp,
                        data = jp2_replication) 

# Generate predictions
predictions_women <- ggpredict(model_women_concern, terms = c("conditions", "female"))
predictions_women$outcome <- "Concerns with Immigrant Women"

predictions_men <- ggpredict(model_men_concern, terms = c("conditions", "female"))
predictions_men$outcome <- "Concerns with Immigrant Men"

predictions_all <- rbind(predictions_women, predictions_men)

predictions_all$outcome <- factor(predictions_all$outcome,
                                  levels = c("Concerns with Immigrant Women",
                                             "Concerns with Immigrant Men"))
predictions_all$x <- factor(predictions_all$x,
                            levels = c("Control",
                                       "Treat1. Skilled Immigration",
                                       "Treat2. Unskilled Immigration"))

# Plot
figure6 <- ggplot(predictions_all,
                aes(x = outcome, y = predicted,
                    color = group, shape = group, linetype = outcome,
                    group = interaction(outcome, group))) +
  geom_point(position = position_dodge(width = 0.6), size = 5, na.rm = TRUE) +
  geom_linerange(aes(ymin = conf.low, ymax = conf.high),
                 position = position_dodge(width = 0.6), linewidth = 1.1, na.rm = TRUE) +
  labs(y = "Predicted Values of Competition Concerns",
       x = NULL,
       color = "Participant's Gender",
       shape = "Participant's Gender",
       linetype = "outcome") +
  scale_color_manual(values = c("#FF8C00", "#800000"),
                     labels = c("Japanese Men", "Japanese Women")) +
  scale_shape_manual(values = c(17, 15),
                     labels = c("Japanese Men", "Japanese Women")) +
  scale_linetype_manual(values = c("solid", "dotted"),
                        labels = c("Concerns with Immigrant Women", "Concerns with Immigrant Men")) +
  facet_wrap(~x, scales = "free_x", nrow = 1) +
  theme_bw() +
  theme(
    strip.text = element_text(size = 13, face = "bold"),   
    axis.title.y = element_text(size = 13, face = "bold"),
    axis.text.y  = element_text(size = 13, face = "bold"),
    axis.title.x = element_text(size = 13, face = "bold"),
    axis.text.x  = element_blank(),
    axis.ticks.x = element_blank(),
    legend.position = "bottom",
    legend.box = "vertical",
    legend.text = element_text(size = 13, face = "bold"),
    legend.title = element_text(size = 13, face = "bold")
  ) +
  guides(
    linetype = guide_legend(title = "outcome", order = 1,
                            override.aes = list(shape = NA, size = 1)),
    color    = guide_legend(title = "Participant's Gender", order = 2),
    shape    = guide_legend(title = "Participant's Gender", order = 2)
  )

print(figure6)

# ============================================================================
# FIGURE 7: Hiring Preferences (MANUSCRIPT)
# ============================================================================

jp2_combined <- jp2_replication %>%
  mutate(gender_label = ifelse(female == 1, "(b) Female Respondents", "(a) Male Respondents")) %>%
  pivot_longer(cols = c("immigrant.women", "immigrant.men", "native.women", "native.men"),
               names_to = "a11b_variable", 
               values_to = "a11b_value")

jp2_summary_combined <- jp2_combined %>%
  dplyr::group_by(conditions, a11b_variable, gender_label) %>%
  dplyr::summarize(mean_value = mean(a11b_value, na.rm = TRUE), .groups = 'drop')

figure7 <- ggplot(jp2_summary_combined, aes(x = a11b_variable, y = mean_value, fill = conditions)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_text(aes(label = round(mean_value, 1)), 
            position = position_dodge(width = 0.9), 
            vjust = -0.5, size = 4, fontface="bold") +  
  facet_wrap(~gender_label, ncol = 2) +
  labs(x = "Preferred Employees", y = "Mean Number of Positions Allocated") +
  scale_x_discrete(labels = c("immigrant.men" = "Immigrant\nMen", 
                              "immigrant.women" = "Immigrant\nWomen", 
                              "native.men" = "Japanese\nMen", 
                              "native.women" = "Japanese\nWomen")) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, size = 14, face="bold"),
    axis.text.y = element_text(size = 14, face="bold"),
    axis.title = element_text(size = 14, face="bold"),
    strip.text = element_text(size = 14, face="bold"),
    legend.text = element_text(size = 14, face="bold"),
    # Add panel borders
    panel.border = element_rect(color = "black", fill = NA, size = 1),
    strip.background = element_rect(color = "black", fill = "lightgray"),
    # Move legend to bottom
    legend.position = "bottom"
  )

print(figure7)


# ============================================================================
# APPENDIX FIGURES H1 & H2: POST-TREATMENT ANALYSIS
# ============================================================================
# Create subsets based on treatment conditions
jp2.t1.sk <- jp2_replication %>% 
  filter(conditions == "Treat1. Skilled Immigration")

jp2.t2.unsk <- jp2_replication %>% 
  filter(conditions == "Treat2. Unskilled Immigration")

# Appendix H1: Gender of Immigrant Thought of
prop_skil <- 100 * prop.table(table(jp2.t1.sk$imagine.sex))
prop_unsk <- 100 * prop.table(table(jp2.t2.unsk$imagine.sex))

# Create dataframe with correct structure
df.img <- data.frame(
  category = c("Female", "Male", "Other", "Female", "Male", "Other"),
  proportion = c(prop_skil["Female"], prop_skil["Male"], prop_skil["Other"], 
                 prop_unsk["Female"], prop_unsk["Male"], prop_unsk["Other"]),
  dataset = c(rep("Treat1.Skilled Immigrants", 3), rep("Treat2.Unskilled Immigrants", 3))
)

# Set correct colors to match the plot (coral/salmon and teal)
fill_colors <- c("Male" = "#ff7f50",     
                 "Female" = "#66c2a5",    
                 "Other" = "#8da0cb")    

# Set factor levels to match the plot order (Male first, then Female, then Other)
df.img$category <- factor(df.img$category, levels = c("Male", "Female", "Other"))

# Create percentage labels
df.img$proportion_label <- sprintf("%.1f%%", df.img$proportion)

# Create the plot
library(ggplot2)

figure_h1 <- ggplot(df.img, aes(x = dataset, y = proportion, fill = category)) +
  geom_bar(stat = "identity", position = "dodge", width = 0.7) +
  geom_text(aes(label = proportion_label),
            position = position_dodge(width = 0.7), vjust = -0.5, size = 3.5) +
  labs(
    x = "Experimental Conditions",
    y = "Percentage",
    fill = "Immigrant's Gender Thought"
  ) +
  scale_fill_manual(values = fill_colors) +
  scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 25)) +  # Set y-axis scale
  theme_minimal() +
  theme(
    axis.title.y = element_text(size = 13),
    axis.text.x = element_text(size = 11, face = "bold"), 
    axis.title.x = element_text(size = 13), 
    legend.text = element_text(size = 13),
    legend.title = element_text(size = 13)
  )

print(figure_h1)


# Appendix H2: Origin of Immigrant Thought of

# Calculate proportions
prop_table_t1 <- 100 * prop.table(table(jp2.t1.sk$imagine.ori))
prop_table_t2 <- 100 * prop.table(table(jp2.t2.unsk$imagine.ori))

# Create dataframes for each dataset
df_t1 <- data.frame(
  origins_thought = names(prop_table_t1),
  proportion = as.numeric(prop_table_t1),
  dataset = "Treat1.Skilled Immigrants"
)

df_t2 <- data.frame(
  origins_thought = names(prop_table_t2),
  proportion = as.numeric(prop_table_t2),
  dataset = "Treat2.Unskilled Immigrants"
)

# Combine the dataframes
df_combined <- rbind(df_t1, df_t2)

# Define the color map to match the plot
color_map <- c(
  "Southeast Asia" = "#66c2a5",    # Teal/green
  "East Asia" = "#fc8d62",         # Orange/coral
  "Europe" = "#8da0cb",            # Light blue
  "Latin America" = "#e78ac3",     # Pink/magenta
  "Middle East" = "#a6d854",       # Light green
  "Other" = "#bebada",             # Light purple
  "South Asia" = "#ffd92f",        # Yellow
  "North America" = "#e5c494",     # Beige/tan
  "Africa" = "#b3b3b3"             # Gray
)

# Set factor levels in reverse order so largest values appear at top
df_combined$origins_thought <- factor(df_combined$origins_thought, 
                                      levels = c("Africa", "North America", "South Asia", 
                                                 "Other", "Middle East", "Latin America",
                                                 "Europe", "East Asia", "Southeast Asia"))

# Format proportions to one decimal place
df_combined$proportion_label <- sprintf("%.1f%%", df_combined$proportion)

# Load required libraries
library(ggplot2)
library(forcats)

# Create the plot
figure_h2 <- ggplot(df_combined, aes(x = proportion, y = origins_thought, fill = origins_thought)) +
  geom_bar(stat = "identity", width = 0.6) +
  geom_text(aes(label = proportion_label),
            hjust = -0.1, size = 3, fontface = "bold") +
  facet_wrap(~dataset, scales = "free_x") +
  labs(
    x = "Percentage",
    y = "Immigrant's Origin Thought"
  ) +
  scale_fill_manual(values = color_map) +
  scale_x_continuous(limits = c(0, 80), breaks = seq(0, 60, 20)) +  
  theme_minimal() +
  theme(
    axis.title.y = element_text(size = 13),
    axis.text.y = element_text(size = 13), 
    axis.text.x = element_text(size = 13), 
    axis.title.x = element_text(size = 13),
    legend.position = "none",  
    strip.text = element_text(size = 13, face = "bold"),
    plot.margin = margin(t = 1, r = 0.5, b = 1, l = 0.5, unit = "cm"),
    panel.grid.major.y = element_blank(),  
    panel.grid.minor = element_blank()     
  )

print(figure_h2)


# ============================================================================
# APPENDIX I: REGRESSION RESULTS TABLES (I1 & I2)
# ============================================================================

library(stargazer)

# Appendix I1: Results for Concerns over Competition with Immigrant Women (OLS)
model_women_concern_base <- lm(w_cncrn_r ~ conditions + female +
                                 inc_11_r +
                                 edu_8_r +
                                 age_5_r +
                                 pride_r +
                                 econ_nation_notsat_r + econ_self_notsat_r +
                                 job_conf_7_r + 
                                 pid_ldp,
                               data = jp2_replication) 

stargazer(model_women_concern_base, model_women_concern, 
          title="",
          header=FALSE, 
          single.row = TRUE, 
          no.space = TRUE, 
          font.size = "small",
          dep.var.labels = "Concenrs over Job Competition with Immigrant Women",
          order = c("conditionsTreat2. Unskilled Immigration:female",
                    "conditionsTreat1. Skilled Immigration:female", 
                    "conditionsTreat2. Unskilled Immigration",
                    "conditionsTreat1. Skilled Immigration",
                    "female",
                    "inc_11_r",
                    "edu_8_r", 
                    "age_5_r",
                    "pride_r",
                    "econ_nation_notsat_r",
                    "econ_self_notsat_r",
                    "job_conf_7_r",
                    "pid_ldp"),
          covariate.labels = c("Treat2(Unskilled):female",
                               "Treat1(Skilled):female",
                               "Treat2(Unskilled.Immi)",
                               "Treat1(Skilled.Immi)",
                               "Female",
                               "Income", 
                               "Education",
                               "Age",
                               "Nationality.Pride",
                               "National.Economy.Dissatisfied",
                               "Own.Finance.Dissatisfied",
                               "Job.Security",
                               "pid.LDP"),
          flip = TRUE)


# Appendix I2: Results for Concerns over Competition with Immigrant Men (OLS)
model_men_concern_base <- lm(m_cncrn_r ~ conditions + female +
                                 inc_11_r +
                                 edu_8_r +
                                 age_5_r +
                                 pride_r +
                                 econ_nation_notsat_r + econ_self_notsat_r +
                                 job_conf_7_r + 
                                 pid_ldp,
                               data = jp2_replication) 

stargazer(model_men_concern_base, model_men_concern, 
          title="",
          header=FALSE, 
          single.row = TRUE, 
          no.space = TRUE, 
          font.size = "small",
          dep.var.labels = "Concenrs over Job Competition with Immigrant Men",
          order = c("conditionsTreat2. Unskilled Immigration:female",
                    "conditionsTreat1. Skilled Immigration:female", 
                    "conditionsTreat2. Unskilled Immigration",
                    "conditionsTreat1. Skilled Immigration",
                    "female",
                    "inc_11_r",
                    "edu_8_r", 
                    "age_5_r",
                    "pride_r",
                    "econ_nation_notsat_r",
                    "econ_self_notsat_r",
                    "job_conf_7_r",
                    "pid_ldp"),
          covariate.labels = c("Treat2(Unskilled):female",
                               "Treat1(Skilled):female",
                               "Treat2(Unskilled.Immi)",
                               "Treat1(Skilled.Immi)",
                               "Female",
                               "Income", 
                               "Education",
                               "Age",
                               "Nationality.Pride",
                               "National.Economy.Dissatisfied",
                               "Own.Finance.Dissatisfied",
                               "Job.Security",
                               "pid.LDP"),
          flip = TRUE)
